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We present large scale 3D particle-in-cell (PIC) simulations to examine particle 
energization in magnetic reconnection of relativistic electron-positron (pair) plasmas. 



The initial configuration is set up as a relativistic Harris equilibrium without a guide 
field. These simulations are large enough to accommodate a sufficient number of 



tearing and kink modes. Contrary to the non-relativistic limit, the linear tearing in- 
stability is faster than the linear kink instability, at least in our specific parameters. 
We find that the magnetic energy dissipation is first facilitated by the tearing insta- 
bility and followed by the secondary kink instability. Particles are mostly energized 
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trie fields produced by the outflows from reconnection. Secondary kink instability 
leads to additional particle acceleration. Accelerated particles are, however, observed 
to be thermalized quickly. The large amplitude of the vertical magnetic field resulting 
from the tearing modes by the secondary kink modes further help thermalizing the 
non-thermal particles generated from the secondary kink instability. Implications of 
these results for astrophysics are briefly discussed. 
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I. INTRODUCTION 

Magnetic reconnection in collisionless relativistic electron-position (pair) plasma plays 
an essential role in a number of important astrophysical problems, such as relativistic jets 
from supermassive black holes, relativistic winds from pulsars and possibly outflows from 
Gamma-ray Bursts. In the Crab pulsar wind, relativistic reconnection in pair plasmas has 
been considered as a primary mechanism to convert magnetic energy to particle energy. [1, 2] 
Although magnetic reconnection has been extensively studied, the physics of fast reconnec- 
tion and particle energization in 3D relativistic reconnection is still a mystery. Magnetic 
reconnection of pair plasmas provides a unique opportunity to examine critically the funda- 
mental reconnection physics due to the absence of Hall term, [3] which is generally thought 
to be the main mechanism for fast reconnection in electron-ion plasmas. Understanding of 
relativistic magnetic reconnection of pair plasmas will also provide useful physical insights 
for both new laboratory experiments of pair plasmas [4] and observations of the high energy 
Universe. 

Previous reconnection studies of relativistic pair plasmas have focused on the fast recon- 
nection and particle acceleration mostly in 2D [1, 5-7] or reconnection onset in 3D with 
relatively small system size. [8-10] It has been shown that instabilities in a cross-field plane 
such as the Kelvin- Helmholtz instability (KHI) [11, 12] and the drift kink instability (DKI) 
[1, 13] are of critical importance to understanding magnetic energy dissipation. [1] The past 
studies used small-scale particle-in-cell (PIC) simulations, (e.g., the length in the outflow 
direction is ~ 25 inertial lengths), however, the characteristic system size in most astrophys- 
ical/laboratory systems is vastly larger. Thus, it is important to study these processes in a 
truly 3D configuration. Fast reconnection in the relativistic limit in 2D has been reported 
to be mostly facilitated by pair pressure tensors. [6, 14] Yin et al. [15] performed the largest 
3D PIC simulations to date of non-relativistic pair plasma reconnection and found that elon- 
gated current sheet is unstable to secondary kink instability and leads to the bending of the 
current sheet. For particle acceleration, past 2D simulations of reconnection mostly started 
from an initially externally imposed reconnection electrical field and showed that particles 
are primarily accelerated near the X lines by electrical fields. [5, 7, 16, 17] However, electri- 
cal fields within a small fraction of the total volume at X lines are difficult to account for 
the large number of accelerated electrons of the total volume. [18] Recently Chen et al. [18] 



reported some surprising satellite observation evidences that energetic electron flux peaks 
at sites of compressed density within islands during electon-ion reconnection in the Earth's 
magnetosphere. Hesse and Zenitani [14] reported that for relativistic pair plasmas in 2D 
without a guide field, most of particle acceleration happens inside magnetic islands possibly 
like a betatron (but no details are provided) (also see Dmitruk, Matthaeus and Seenu [19]). 
Thus, the well accepted particle acceleration mechanism needs to be revisited, especially in 
3D. Drake et al. [20] suggests that Fermi acceleration from contracting magnetic islands is 
the major mechanism for acceleration in 3D non-relativistic electron-ion plasma magnetic 
reconnection and electrons are mostly accelerated within the islands. In this paper, we will 
present the first numerical study of 3D magnetic reconnection of relativistic pair plasma 
with an emphasis on particle energization. 

The paper is organized as follows. Section II describes the computational model and 
problem setup. In Sec. Ill, overall simulation evolution is briefly described followed by par- 
ticle energization and its mechanisms, particle thermalization and the role of the secondary 
kink instability on particle acceleration in magnetic reconnection. Conclusions and their 
implications for astrophysical situations are given in Sec. IV. 

II. SIMULATION SET-UP 

The geometry of the simulation is given in Fig. 1. The simulations are initialized with 

a relativistic Harris equilibrium current sheet. [2, 7] The initial magnetic field has B y = 

B tanh(x/L) and B x = B z = 0. Here L is the half thickness of the current sheet. The 

initial plasmas consist of two distinct parts: 

, ( ,_ n cosh~ 2 (x/L) 7a(g- (3 s m s cu z ) n^ e_ 
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(1) 
where the subscript s denotes positron (i) and electron e, respectively. u z is the ^-component 

of the four-speed, e is the energy of the particle and K 2 (x) is the modified Bessel function 

of the second kind. One part (the underlined term) is a uniform background plasmas with 

density rib = 0.3no and temperature T\, = T e = Tj = m e c 2 where c is the speed of light and 

m e = rrii = m are the mass of electron and position, respectively. The second part of plasmas 

(for both electrons and positrons) has the same temperature T = m e c 2 as the background 

plasma but their density is given as n = no cosh~ 2 (:r/L). Electrons and positrons are drifting 



against each other, with electrons having +V^ and positrons having — V^. The drift velocity 
Vd is given as /3 = Vj/c = 0.82. Both species have relativistic Maxwellian distributions. 
The initial current from the drifting distributions is in the z direction. The simulation 



parameters are L/di = 0.7 and oj p /Vl c = 0.5, where di = c/u> p , co p = ^ {AnriQe 2 ) / rrii and 
Q c = eBo/(rriic) are the positron inertial length, plasma frequency and cyclotron frequency, 
respectively. Periodic boundary conditions are used in y and z with reflecting boundaries in 
x. 

In such anti-parallel geometry, the linear Vlasov theory predicts two types of instabilities: 
tearing, with wave vectors along y, and kink, along z. [23] Using the parameters given above, 
linear theory predicts that the fastest growing tearing mode has wave number k y L = 0.4 
and growth rate T/Q c = 0.069, while the fastest growing kink mode has k z L = 0.25 and 
r/Q c = 0.059. We have also performed 2D relativistic nonlinear tearing simulation (one 
cell in z direction) and 2D relativistic nonlinear kink simulation (one cell in y direction). 
We find that the fastest tearing mode has k y L = 0.4 and T/Q c = 0.063, and the fastest 
kink mode has k z L = 0.25 and T/Q c = 0.053. They match the linear Vlasov theory very 
well. Petri and Kirk [21] also carried out similar linear Vlasov analysis in the relativistic 
limit (see also [22]). Compared to the results in the non- relativistic regime, [23] the growth 
rates of both modes are greatly reduced and the wavelengths of both modes are increased 
by the relativistic effects, which is further verified in the 3D simulations. While in the non- 
relativistic regime the growth rate of the fastest growing kink mode is almost twice as large 
as that of the fastest growing tearing mode, the growth rate of the fastest tearing modes in 
the relativistic regime is a bit larger than the fastest growing kink modes, which results in 
the earlier emergence of tearing modes than kink modes in 3D simulations, contrary to that 
in the non-relativistic regime. 

In 3D, to examine how the reconnection dynamics may be influenced by both the lin- 
ear and nonlinear competition between the tearing instability (in x-y plane) and the kink 
instability (in x-z plane), we use four simulations with the same size in the tearing plane 
L x x L y = 200di x 200di but different sizes in the z direction so that different spectra of kink 
modes can be excited. They include L z = 0.195rfj (only 2D tearing is present), 20di, 50d i} 
and 200dj, which are referred to as run A, B, C and D, respectively. These simulations follow 
the dynamics of 0.168, 16, 43, and 84 billion particles on 1024 x 1024 x 1, 1024 x 1024 x 96, 
1024 x 1024 x 256, and 1024 x 1024 x 500 meshs, respectively. For these simulations, the cell 



sizes are: 5x = 5y = 0.195Ad, while Sz is 0.195A/), 0.208Ad, 0.195A/), and 0A\ D , respec- 
tively. Here, \ D is the Debye length, and the time step 5tVt c is 0.223, 0.228, 0.223, and 0.258, 
respectively. The total energy is conserved to better than 0.1% throughout the simulation 
runs. In order to remove the recirculation issue, all data presented here are chosen before 
recirculation happens (^^recirculation ~ QcL y /v gw ~ 4,000, where v 9jy ~ 0.1 is the typical 
particle group velocity in y-direction). 

III. RESULTS 

A. Overall Evolution 

We now discuss the dynamics of the evolution. 

Fig. 2 shows the time evolution of total magnetic energy of the four simulations. We 
have divided the magnetic energy conversion evolution into two main stages (separated by 
the dash line in Fig. 2). For the first stage (tearing stage thereafter) when tQ c < 954, the 
dynamics is dominated by the development of relativistic tearing instability as indicated by 
the red and green curves. The linear kink instability plays essentially little role for energy 
conversion in run A and B. It becomes visible in run C and D but it seems subdominant to the 
tearing modes (it has smaller growth rate compared to the tearing mode). The purely linear 
growth stage lasts until tVt c ~ 250, then the islands from tearing start to merge, entering the 
nonlinear stage. The time tfl c ~ 954 marks the eventual saturation of the tearing nonlinear 
evolution. For the second stage (secondary kink stage thereafter) when tQ c > 954, a new 
nonlinear instability becomes active. This was referred to as the secondary kink instability 
in an earlier, similar study using non-relativistic pair plasmas. [15] This happens in 3D 
only when there is sufficient length in the z-direction, as indicated by run C and D (black 
and blue curves). We can see that the tearing instability converts about 27% of the initial 
magnetic energy into particle energy, followed by the additional conversion from ~ 27% to 
~ 40% of the initial magnetic energy by the secondary kink for run C and D. This result 
indicates that the overall energy conversion depends on how many modes are allowed in the 
z-direction and their nonlinear interactions with the tearing instability. 

The relativistic cases presented here show a number of similarities to the previous non- 
relativistic 3D cases studied in Yin et al. [15]. During the linear phase tfi c < 250 (Fig. 3), 



about 14 small magnetic islands form along the (x, y) plane from the tearing instability. 
Contrary to the non-relativistic limit, [15] the primary kink mode with wavelength ~ 18.5rfj 
appears a little bit later than the tearing mode, undulating the small magnetic islands in the 
z direction (Fig. 3). Then the magnetic islands formed from the tearing instability coalesce to 
produce a dominant reconnection site through which most of the magnetic flux is processed 
(tQ c ^ 954) (Fig. 4). When magnetic islands are formed, the islands are not observed to be 
contracting, but rather they are either expanding or merging or remaining in a quasi-steady 
state. For run D, we see that, in the linear phase, both the tearing and kink instabilities 
are growing, causing the the reconnection sites to be "patchy" [Fig. 4 (c)]. After the major 
reconnection site is established (tQ c > 954) (Fig. 5), this dominant diffusion region expands 
in y to form a highly elongated current sheet unstable to secondary kink instability if the 
dimension in the current direction (z) is large enough. Note that the wavelength ~ 40<ij of 
the secondary kink mode is larger than the wavelength of the linear kink mode, thus there 
is no secondary kink modes in run B (L z = 20di). The current sheet is bent by secondary 
kink instability and form plasmoids. In the case of L z = 20di (run B), the elongated current 
sheet is also unstable to the secondary tearing instability, small secondary magnetic islands 
appear, expand and merge with the large primary magnetic islands. In all cases this major 
reconnection site becomes quite extended in the z direction (Fig. 5). 

An analysis of the generalized Ohm's law reveals that the reconnection electrical field is 
mostly from pair pressure tensors at the reconnection sites (Fig. 6), consistent with previous 
results. [3, 14] The dimensionless reconnection rate Er =< E z /B > (<> represents an 
average over the inflow surface, vertical electrical field component E z and magnetic field B 
are measured at the inflow surface) after the establishment of the major reconnection site 
remains fast and time dependent, with Er ps 0.068 and 0.064 for run B and C, respectively. 
These rates are smaller than the value in the non-relativistic regime. [15] The reduction of 
reconnection rate might be due to the higher effective inertia in the relativistic case. [7] The 
inflow surface is chosen to be ±23dj upstream from the dominant reconnection site and has 
an extent of 12.5d{ in y. The time- dependent nature of the reconnection rate is due to both 
the formation and destruction of the plasmoids (similar to the non-relativistic case). 



B. Particle Energization: Time and Location 

We now describe the particle energization processes. Given the fact that the 2D run gives 
a similar amount of total energy conversion as 3D runs during the tearing stage, we first 
discuss the magnetic dissipation process associated mostly with the tearing instability in the 
{x, y} plane and assume that the variation along the z direction does not play a major role. 
During the second stage, for 3D runs with sufficient length in the z direction, the variation 
along the z direction plays an important role in particle energization (see §111 E below). 

Fig. 7 shows the particle energy distributions at different times based on two runs, A 
(panels a, b, and d) and C (panel c). Panels a and c represent the particle distribution 
from the whole computational domain of run A and C. Panels b and d represent the particle 
distribution taken from a sampling box with Ax = Ay = Id^ at the magnetic island (b) and 
the reconnection site (d) of run A. The four curves in each panel represent four characteristic 
stages of the evolution: the initial state, the end of linear phase, the end of nonlinear 
interaction stage, and during the secondary kink (for run C), respectively. Note that, for 
panels b and d, the total number of particles within the sampling box is changing with time. 
We have normalized the distributions at different times so that the curves have the same 
number of particles within each panel. 

The initial particle distributions [the black curves in panel (a, c)] are composed of two 
parts: one part from the background plasmas and the other from the current sheet plasmas, 
which has a higher effective temperature because of the drift. On the other hand, the black 
curves in panel (b, d) are mostly composed of the current sheet plasmas only. Over all, a 
major part of energy conversion takes place between the end of the linear phase and the 
establishment of the major reconnection site (between blue and green lines in panels a and 
c), consistent with Fig. 2. 

Inspecting panel d, we find that there is some amount of particle acceleration at the 
reconnection site, as indicated by the blue curve. This is consistent with the fact that the 
reconnection electric field is the largest at these locations. But such acceleration, although 
having large local electric fields, is very limited spatially and does not last in time since the 
particles can not stay at the reconnection site for a long time but rather quickly leave the 
reconnection sites via outflows and enter the magnetic islands, as indicated by the green 
and red curves when the particle energy distribution becomes less energetic than the initial 
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distribution. Eventually the plasmas at the reconnection sites are replaced by the edge 
plasmas from inflow, from the black line with nonzero (3 to the red line with zero (3. (Note 
that for a given drifting relativistic Maxwellian distribution, the absolute value of the slope 
is proportional to oc (1 — /3)/T using the log-linear scale, where /3 is the drift velocity.) 

Panel b, however, reveals that a large amount of energy conversion occurs in the magnetic 
island region. Even though the electric field in these island regions might not be as high when 
compared to the peak reconnection electric field, both the volume and number of particles 
that experience these electric fields are much larger. The particles are confined inside the 
magnetic islands for a much longer time than at the reconnection site. Furthermore, by 
analyzing particle distribution changes within much smaller time spans (see below), we find 
that the energetic particle population in the island regions must be produced in-situ, i.e., 
the majority of the energetic particles do not get energized in the reconnection region then 
propagate to the island regions. Note that plasmas inside the islands seem to be eventually 
thermalized into another relativistic Maxwellian distribution with a temperature roughly 2.3 
times the initial temperature of the background plasmas [the green and red lines of Fig. 7 
(a,b)]. 

C. Particle Energization: Mechanism 

Although particle acceleration is expected at the reconnection sites, it is not so obvious 
why the majority of particle energization occurs in the island regions in these simulations. 
During the tearing stage, outflows from reconnection (cf. Fig. 1) lead to strong and y- 
dependent E — — (VxB)/caX the edges of islands (see Fig. 1), where V is mostly along the 
y-direction and B is in the {x, y} plane. The resulting electric field component E z is quite 
significant because the outflow velocity can become quasi-relativistic. This is a peculiarity of 
the relativistic pair plasmas because of the absence of ions, which are much heavier than the 
electrons. Therefore the Alfven velocity (the typical outflow speed for particles) can become 
relativistic too, leading to a large component of —V x B. [8] We find that these electric fields 
tend to be long-lived and are distributed over large spatial volumes. Consequently they 
are responsible for the majority of magnetic energy conversion to particles, at least for the 
configurations we have studied here. 

To understand the details of particle energization, we plot the spatial distribution of the 
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normalized E 2 — B 2 in Fig. 8, using results from run A. The reconnection regions indeed have 
E 2 — B 2 > 0, indicating the existence of a net electric field that can accelerate particles. 
Most other regions have E 2 — B 2 < 0, which has been usually interpreted as not useful 
for particle acceleration because there exists a Lorentz transformation that makes electric 
field vanish. [14] Such consideration, however, does not apply when there are large spatial 
variations of E z or B x or both. In the bottom panel of Fig. 8, we plot the B x and E z 
profiles along the |/-axis (with x — 0). As particles try to gyrate in the {y, z} plane due to 
B x but experience a spatially varying E z along the ^-direction within their gyro-orbits, net 
acceleration can occur. 

Particular types of field structures are favored in particle acceleration as shown in Fig. 9. 
Panel (a) emphasizes the spatial variations in E z . Particles gain more energy in the larger 
\E Z \ side than the energy lost in the smaller \E Z \ side. Another scenario emphasizes the 
variations in B x . Particles have larger gyroradius in the smaller \B X \ side than in the larger 
\B X \ side, leading to net particle acceleration in the z direction. [16] Fig. 8 {bottom) shows 
that both of these structures are present when the tearing instability becomes nonlinear. 
To confirm this acceleration process, we have performed test particle calculations using the 
initial particle distribution as in the simulation [2, 7] and with three different field structures: 
1) constant B x but spatially varying E z [Fig. 9 (a)]; 2) constant E z but spatially varying 
B x ; 3) temporally and spatially varying B x and E z directly from the 2D simulation run 
A [Fig. 9 (&)]. The E z and B y profile data from the nonlinear simulation are fitted with 
first order polynomials by the least square fit method at tfl c ~ 329 between y = 60di and 
y = HOdi. And then a spatially constant E z — 0.1 and a realistic B y profile are used for 
(1) and a spatially constant B y = 1.0 and a realistic E z profile are used for (2). The total 
number of particles used in these test particle calculations is around 500, 000. We find that 
the varied E z scheme [Fig. 9 (a)] produces ~ 10 times more efficient particle acceleration 
than the varied B x scheme. Quantitatively, for efficient acceleration, it is necessary for the 
typical E z variation scale to be smaller than the particle's gyroradius: 

\E z /(dE z /dy)\<^mc 2 /(e\B x \), (2) 

where E z is the characteristic E z along the particle trajectories. This suggests that the 
more sharply varied E z or weaker B x leads to more efficient particle acceleration, which is 
confirmed by our test particle calculations [Fig. 9 (c)]. Furthermore, test particle calculation 
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using the time-dependent field output from the simulation run A gives d^y/d(tQ c ) ~ 0.02 by 
statistical average [Fig. 9 (d)], which is sufficient to explain the particle energy gain before 
the secondary kink instability kicks in. We confirm that particles are mostly accelerated at 
the edges of the islands (e.g., the region around y ~ 80 and 90rfj, also see Fig. 1). This type of 
field configuration is quite robust and is found to persist throughout the whole tearing stage. 
The long-term (tfl c > 500) test particle calculations result in quasi-Maxwellian distributions 
at low energy, but with significant (> 15%) non-thermal high-energy particles (7 > 20) at 
tVt c > 750 [Fig. 9 (d)], which is not observed in the real simulations. This is most likely due 
to the thermalization process in real simulations (see below). 

D. Particle Energization: Thermalization 

The test particle study shown in Fig. 9 suggests that the fields from the tearing instability 
could accelerate particles to high energies (7 ~ 100). Fig. 10, however, indicates that 
the non-thermal high-energy particles from this acceleration process have been significantly 
thermalized. The particle distributions are taken from a sampling box at x = and y = 89rfj, 
which is near the right edge of the central island shown in Fig. 8. Specifically, as electrons 
are accelerated in the +p z direction (panel b of Fig. 10), their +p z momentum turns into 
— p y momentum due to the relatively large B x , e.g., the extended tails in — p y shown as 
unsymmetric black dash line in panel a. Note that the positrons will experience similar 
processes, except that they have — p z turning into —p y . Subsequent evolution shows that 
there exists a series of plateaus with decaying magnitudes at the — p y side of particle number 
distribution [blue and green lines of Fig. 10 (a)]. The particle distribution then turns into an 
approximately symmetric Maxiwellian distribution in p y as shown by the red line of Fig. 10 
(a). This distribution can be fitted with two temperatures, suggesting that the energy gained 
through the E z acceleration has been thermalized in the p y direction. 

Note that the time difference between the black dash and red curves in Fig. 10 (a) is 
5tQ c ~ 50, this implies that the thermalization process is quite efficient. We suggest that a 
mixture of bump-on-tail and two-stream instabilities can cause wave excitation and Landau 
damping, and this eventually results in a symmetric distribution in p y . This process could 
be fast because the strength of two-stream instability and Landau-damping is proportional 
to plasma frequency u p in pair plasmas , much larger than in electron-ion plasmas (~ 
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ujp e ^m e /mi). The above particle thermalization process would happen repeatedly inside 
the islands and eventually all the plasmas inside the islands are thermalized into a higher 
temperature. From Fig. 10 (a), the typical thermalization time scale is found to be StQ c ~ 
50. 

It is worth noting that during the time window shown in Fig. 10, no island merging 
is observed to take place. Furthermore, similar particle distribution evolution shown in 
Fig. 10 appears also at both downstream (y = 87dj) and upstream (y = 91di) regions 
simultaneously around y = 89di region. This rules out the possibility that this energization 
is due to the particle convection, but supports the suggestion that the thermalization is due 
to local "micro-instabilities" . In addition, we observe that particles are drifting slowly in the 
y direction, with a typical particle velocity v y /c = Py/i'ym) < 0.2 (not phase velocity). This 
means that the particles travel < 5di in StQ c ~ 50 in the ^/-direction, which is much smaller 
than the size of the magnetic islands (~ 30<ij), thus particles do not have enough time to 
bounce between the edges of the islands. The typical gyro-radii of the particles are always 
5 — 7 times smaller than the island size at every stage, thus the particle acceleration from 
bouncing between "walls" of contracting magnetic islands [20] due to the fast gyro-motion 
is not possible either. 

We also find that particles have very small p x throughout the evolution and the p z distri- 
bution remains "skewed" due to the presence of E z until quite late time (tQ c ~ 1800) when 
the distribution in p z becomes symmetric as well [Fig. 10 (&)]. 

Even though we mostly present the results from run A, we found similar results with 
run C. Both the acceleration mechanism that relies on the spatial variation of E z and the 
subsequent thermalization process seem to be robust processes for particle energization in 
pair plasmas. 

E. Particle Energization: The Role of the Secondary Kink Instability 

For both run C and D, there is additional particle energization in the secondary kink stage. 
We find that current sheets along the z direction emanating from both the reconnection 
and island regions undergo the secondary kink instability with similar wavelength ~ 40rfj 
(Fig. 11). The magnitude of kinking at the island region is larger because the kinking 
at the reconnection site is suppressed by the reconnection inflow (Fig. 11). The onset of 
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the secondary kink instability inside the island is also delayed compared to the one at the 
reconnection site since the the current sheet at the magnetic islands is thicker. Those 
secondary kink instabilities are new kink modes since the the initial primary kink modes 
have been stabilized during the tearing stage. 

The appearance of the z-dependent B y and E z resulting from the secondary kink insta- 
bility (Fig. 11) looks very similar to Fig. 1 of Zenitani and Hoshino [1], which suggests that 
the same particle acceleration mechanism works for the secondary kink instability during 
this stage, i.e., a "DC" acceleration channel mechanism (for detailed discussion please see 
Zenitani and Hoshino [1]). This same mechanism accelerates high-energy particles (7 > 10) 
efficiently at both magnetic islands and reconnection sites. This leads to the production 
of the non-thermal particles [high energy tail of red line of Fig. 7 (c)], although those 
non-thermal particles would be mostly thermalized by the above thermalization process as- 
sociated with B x and only some of non-thermal particles might be left. Here E z ~ 0.2 is 
mostly due to the charge separation from the relative electron and positron drift, [7] which 
leads to the onset of the secondary kink instability. 

The high-energy particles undergo a fast streaming motion in the z direction but at the 
same time executing frequent cyclotron motion in the (p y ,p z ) space, leading to efficient ther- 
malization as discussed in Sec. Ill D. Such a strong thermalization process is not observed 
in previous 2D kink studies [1] since B x ~ 1 induced from tearing modes in this study is 
much larger than in those cases ~ 0.1. This might explain a relatively higher percentage of 
non-thermal particles found in those studies [1]. 

It is worth emphasizing that, at least for the initial configuration we have studied here, 
the particle energization by tearing occurs mostly inside the magnetic islands and produces 
mostly thermal heating in both stages. But particle energization by the secondary kink 
instability leads to some non-thermal particles in the second kink stage. These effects 
could only be studied with 3D simulation of a large enough computational domain in every 
direction. 

IV. SUMMARY 

In this paper we have presented large scale 3D particle-in-cell (PIC) simulations to ex- 
amine particle energization in magnetic reconnection of relativistic electron-positron (pair) 
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plasmas. These simulations are large enough to accommodate a sufficient number of tearing 
and kink modes. The interplay of these two main instabilities regulates the overall evolution: 
from the initial linear tearing and kink growth with the appearance of small magnetic islands. 
to the merging of small magnetic islands into a big one, establishing the major reconnection 
sites, to the elongated current sheet becoming unstable to the secondary kink instability, 
forming plasmoids. The particle energization could be divided into two main stages: (1) 
the tearing stage (tQ c < 954), and (2) the secondary kink stage (tQ c > 954). We find that 
particles are mostly energized inside the magnetic islands during the tearing stage due to 
the spatially varying electric fields produced by the outflows from reconnection. Secondary 
kink instability also leads to some particle acceleration. Accelerated particles are, however, 
observed to be thermalized. Since the physical conditions for both the energization process 
and the thermalization process discussed in this study are quite typical in relativistic pair 
reconnection, we suggest that these processes are germane in understanding the magnetic 
energy dissipation in pair plasmas. 

In a relativistic current sheet, previous literature reported that the linear kink instability 
is faster than the linear tearing instability unless the reconnection is driven. [1, 7, 9, 10] In this 
work, contrary to the above well-known results, the linear tearing instability is faster than 
the linear kink instability, at least in our specific parameters. This result raises another hint 
to the striped wind problem. [2] It is worth of emphasizing that all the calculations reported 
in this paper are done with an initial anti-parallel magnetic field without a guide field for pair 
plasmas. The introduction of a uniform guide field would significantly change the results as 
reported in Hesse and Zenitani [14]. Compared to a relativistic reconnection without a guide 
field, the reconnection rate is slower with a guide field due to a lower compressibility effect 
of the guide magnetic field. More importantly with a guide field, the particle acceleration 
is mostly done at the reconnection site, whereas it is inside the magnetic island without 
a guide field. [14] The magnetic field configuration in the real astrophysical situation such 
as pulsar wind would be even more complicated. We need to be cautious when trying to 
apply the particle energization processes reported in this paper to the real astrophysical 
situations. Additionally the energization process in pair plasmas may find only limited 
applications in the solar or magnetospheric environments where electron-ion reconnection 
should be dominant. 
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Figure 1: (Color) Overall geometry of the simulation. Two main types of instabilities are excited 
in this configuration: the tearing instability along the k y direction and the kink instability along 
the k z direction. Initially there is an anti-parallel magnetic field (cyan arrows) in the x-y plane 
and a current sheet (pink) in the z direction at the interface. After an initial perturbation from 
thermal noise, magnetic reconnection in x-y plane would take place. Inflow (blue vertical arrows), 
outflow (blue horizontal arrows) and magnetic islands would form. Different from 2D magnetic 
reconnection, in the z direction, the strong current induces current-driven instabilities such as kink 
instability. The edges of magnetic islands are also indicated by black arrows. 
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Figure 2: (color) Time evolution of magnetic energy with L x = L y = 200dj. Red, green, black 
and blue curves are for L z = 0.195 (2D), L z = 20,50, and 200dj, respectively. All energies are 
normalized to the initial total magnetic energy. The current sheet starts to bend (kink unstable) 
in the z direction at tft c ~ 954 (dash line). 
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Figure 3: (Color) Linear phase of the simulation. Many small magnetic islands resulting from 
tearing undulate in the z direction due to primary kink instability (tD, c < 250). The blue-red 
regions are iso-surfaces of density with color indicating peak of B y (blue: negative; red: positive). 
These iso-surfaces enclose the volume of the magnetic islands. Embedded are the green regions 
for the iso-surfaces of weak \B\ with color indicating peak of \J\. These iso-surfaces enclose the 
volume of the current sheets at the reconnection sites. Contrary to the non-relativistic limit, the 
tearing mode appears earlier than the kink modes, (a): run B,L Z = 20^; (6): run C, L z = 50dj; 
(c): runD, L z = 200^. 
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Figure 4: (Color) Nonlinear interaction stage of the simulation. The blue-red regions are iso- 
surfaces of density with color indicating peak of B y (blue: negative; red: positive). These iso- 
surfaces enclose the volume of the magnetic islands. Embedded are the white regions for the 
iso-surfaces of weak \B\ with color indicating peak of \J\. These iso-surfaces enclose the volume 
of the current sheets at the reconnection sites. The magnetic islands formed from the tearing 
instability coalesce to produce a dominant reconnection site through which most of the magnetic 
flux is processed (tO c < 954). The localized patchy reconnection sites evolve as they self-organize 
in z to form a single, large diffusion region. 
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Figure 5: (Color) Saturation stage of the simulation. The blue-red regions are iso-surfaces of 
density with color indicating peak of B y (blue: negative; red: positive). These iso-surfaces enclose 
the volume of the magnetic islands. Embedded are the white regions for the iso-surfaces of weak 
\B\ with color indicating peak of \J\. These iso-surfaces enclose the volume of the current sheets 
at the reconnection sites. The dominant diffusion region expands in y to form a highly elongated 
current sheet unstable to the secondary kink instability. The current sheet is bent by the secondary 
kink instability and form plasmoids (tQ c > 954). 
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Figure 6: (color) x distribution of ~E Z (black solid), [E + V x B/c] 2 (black dash), electron pressure 



term[-^V-P e ], 



(green), positron pressure term [^V • Pj] 2 (yellow), convective bulk inertial 



term [^^^ ' (J v + V J)]^ ( c Y an ) and Jz (blue) at an X point (x,y,z) = (0,232,0) of run B at 
tQ c i ~ 1011. The time dependent term [y^i%;}z, which is hard to calculate accurately, is not shown 
here. [3] All terms are normalized by BqVa/c, where Va = -Bo/[47rno(m e + rrii)}. J z is normalized 
by tiocVa- x is normalized by the half thickness of the initial current sheet L. 
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Figure 7: (color) Particle number distribution v.s. energy at some typical stages of run A and C. 
For run A, panel (a), (b) and (d) represent the particle distributions from the whole computation 
domain, at an magnetic island, and at a reconnection site, respectively. For run C, panel (c) 
represents the particle distribution from the whole computation domain. All curves are normalized 
to have the same total number of particles. For (a,b,d), curves black, t£l c = 0; blue, t£l c = 149; 
green, tQ c = 808; red, tU c = 1811. For (c), curves black, tU c = 0; blue, tQ. c = 238; green, iO c = 954; 
red, tfl c = 1470. 
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Figure 8: (Color) Contour plot of magnetic flux function ^ = — J B y dx colored by E 2 — B 2 . 
The dash contour lines of magnetic flux function ^ defines the magnetic islands structures. Three 
magnetic islands are observed at this moment. Fig. 9 focuses on the middle island. 
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Figure 9: (Color) (a) Schematic view of the acceleration mechanism due to the variation of E z 
(constant B x ). (b) The y profile of B x {black) and E z (red) at t£l c = 329 of run A. (c) Normalized 
particle distribution v.s. energy 7 with different B x strength and a fixed E z profile: blue, B x = 0.7; 
green, B x = 1.0; red, B x = 3.0; black: normalized initial particle distribution for comparison. 
Larger B x leads to higher energy threshold (Eq. 2), therefore fewer low-energy particles are accel- 
erated, (d) Particle distribution evolution v.s. energy 7 using time-dependent field data from the 
simulation run A. black, tQ c = 0; blue, t£l c = 300; green, t£l c = 500; red, t£l c = 750. Significant 
non-thermal high-energy particles are obtained at t£l c ~ 750 (red line). 
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Figure 10: (Color) Particle number distribution v.s. p y (a) and electron number distribution v.s. 
p z (b) of run A, taken from a box with Ax = Ay = 2d% centered at x/di = and y = 89di, at 
the right edge of the central island. All data are normalized, (a): black solid, t£l c = 0; black 
dash, t£l c = 295; blue, t£l c = 311; green, iO c = 329; red, t£l c = 348. (b): black solid, t£l c = 0; 
blue, t£l c = 329; green, ttt c = 425; red, t£l c = 1811. 




Figure 11: (Color) Spatial distribution of plasma density and fields at the nonlinear saturation 
stage from run C with tVt c ~ 954 . There are three cutting planes. One in the {x, y} plane with 
contours of particle density p € [0.09, 0.4] (from pink to white) at a z location near the boundary, 
indicating the elongated reconnection site along y between the large magnetic islands. The other 
two are in the {x,z} plane with contours of reconnection magnetic field B y (white solid lines), 
colored by accelerating electrical field E z S [—0.25,0.24] (from blue to yellow) at magnetic island 
(y/di = 60) and at reconnection site (y/di = 134) respectively. The run D with L z = 200dj gives 
the same conclusions about the second kink instability with similar wavelength. 



